Finite temperature excitations of a trapped Bose gas 
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We present a detailed study of the temperature dependence of the condensate and noncondensate 
density profiles of a Bose-condensed gas in a parabolic trap. These quantitites are calculated self- 
consistently using the Hartree-Fock-Bogoliubov equations within the Popov approximation. Below 
the Bose-Einstein transition the excitation frequencies have a realtively weak temperature depen- 
dence even though the condensate is strongly depleted. As the condensate density goes to zero 
through the transition, the excitation frequencies are strongly affected and approach the frequencies 
of a noninteracting gas in the high temperature limit. 
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Since the recent discovery of Bose-Einstein condensa- 
tion (BEC) in ultracold trapped atomic gases there 
have been a number of calculations of the collective ex- 
citations in such systems. These calculations J3|-[| have 
mainly concentrated on the excitations at T = using 
the standard Bogoliubov equations, which assume that 
all the atoms are in the Bose condensate. Recent ex- 
perimental results for rubidium JtJ and sodium |8| are 
in good agreement with these T = calculations, which 
give the frequencies as a function of the number of atoms. 
However, at finite temperatures where the condensate is 
strongly depleted, a generalization of these theories is 
needed. The simplest theory which includes the effect of 
the noncondensate atoms in a self-consistent manner is 
the Hartree-Fock-Bogoliubov (HFB) approximation B. 

In this Letter, we use a simplified version of the HFB 
approximation originally introduced by Popov [ p^[ in 
the study of uniform Bose gases close to the transition 
and discussed more recently p]Jl as a useful collision- 
less approximation over a wide range of temperatures. 
The Popov approximation ignores the anomalous (or off- 
diagonal) correlations present in the full HFB theory. 
Within this approximation, we calculate self-consistently 
the condensate n c (r) and noncondensate n(r) density 
profiles at T — as well as at finite temperatures. Our 
present results are the first self-consistent calculation 
of these quantities within the HFB-Popov theory which 
avoids the semiclassical approximation for the quasipar- 
ticle excitations used previously jL2j, or various other 
simplifications |0|. Our calculations also provide the 
HFB-Popov excitation frequencies, which at low temper- 
atures are the collective modes of the condensate. We 
find that these modes are weakly temperature-dependent 
in a range of temperatures below T c , even when the con- 
densate is strongly depleted. At temperatures above T c , 
the lowest excitation frequencies approach the excitation 
energies of noninteracting particles in a parabolic well. 

Starting with the equation of motion for the Bose field 
operator ip(r), and its decomposition ip(r) = $(r) + V>(r) 
in terms of condensate and noncondensate parts 



the condensate wavefunction $(r) is defined within the 
Popov approximation by the generalized Gross-Pitaevskii 
equation KM 



2m 



Vext(r) +g[n c (r) + 2n(r)] 



$(r)=/x$(r). (1) 



Here, n c (r) = |$(r)| 2 and n(r) = (ip^ (r)ip(r)) are the 
condensate and noncondensate densities, respectively. 
The Popov approximation [f^llf consists of omitting the 
anomolous correlation (il>(r)ip(r)), but keeping n(r). The 
condensate wavefunction in Eq.(|l|) is normalized to N c , 
the total number of particles in the condensate. V ext (r) 
is the external confining potential and g — 4irh 2 a/m is 
the interaction strength determined by the s-wave scat- 
tering length a. The condensate eigenvalue is given by 
the chemical potential /i |l4j . 

The usual Bogoliubov transformation, ip(r) = 
J2ii u i( r )&i — v i( r )®t]i 1,0 the new Bose operators on and 
a\ leads to the coupled HFB-Popov equations 



Cui(r) - gn c (r)vi(r) 
£vi(r) - gn c (r)ui(r) 



E l u l (r) 
-EiVi(r) , 



(2) 



with L = — V 2 /2m + V ext (r) + 2gn(r) — fi = h + gn c (r). 
These equations define the quasiparticle excitation ener- 
gies Ei and the quasiparticle amplitudes Ui and Vi . Once 
these quantities have been determined, the nonconden- 
sate density is obtained from the expression MM 



n(r) 



ni(r) + n 2 (r) 



(3) 



where hi(r) is that part of the density which reduces 
to the ground-state noncondensate density as T — > 0. 
The component H2(r) depends upon the Bose distribu- 
tion N (Ei) = (e pE ' - l)" 1 and vanishes in the T -> 
limit. 
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Rather than solving the coupled equations in Eq.(g) 
directly, we introduce a new method based on the auxil- 
iary functions ip^'(r) = Uj(r)±«j(r). From Eq.(||), these 
are found to be solutions of the uncoupled equations 

hl4 +) (r) + 2gn c (r)h 4 +) (r) = E?4 +) (r) (4) 

(5) 



h 2 4-\r) + 2gh Q n c (r)^->(r) = E?iPf>(r) . 



(-)/ 



The two functions are related to each other by hoW — 

Eiipi . We note that the collective modes of the con- 
densate can be shown to have an associated density fluc- 
tuation given by Srn(r) oc $(r)?/>^' ) (r). 

To solve Eq.(|) we introduce the normalized eigenfunc- 
tion basis defined as the solutions of ho<f) a (r) = e a <p a (r). 
The lowest energy solution gives the condensate wave- 
function <£>(r) = v / AoV , o( r ) with eigenvalue Eq — 0. Us- 



ing the expansion ip^ + \r) = 7"L Ca'4> a (r) in Eq.(0), we 



obtain the eigenvalue equation |lq] 







{M a0 



(G) 



where the matrix M Q( a = 2g J dr <p^(r)nc(r)(pp(r). The 
orthonormality of the u and v functions |1 in Eq. (j^) re- 
quires that the expansion coefficients c~a be normalised 
as ^ Q e a (| c« = EiSij. The two parts of the non- 
condensate density n(r) as defined in Eq.([s]) are given in 
terms of the Ca through the expressions 



(7) 



and 



xiVo(^). 



(8) 



The calculational procedure can be summarized for an 
arbitrary confining potential as follows: Eq.(|l|) is first 
solved self-consistently for <&(!"), with h{r) set to zero. 
Once $(r) is known, the eigenfunctions of ho required in 
the expansion of the excited state amplitudes are gener- 
ated numerically. The matrix problem in Eq.(^) is then 
set up to obtain the eigenvalues Ei, and the correspond- 

(i) 

ing eigenvectors c a are used to evaluate the nonconden- 
sate density. This result is inserted into Eq.(|l|) and the 
process is iterated to convergence. At each step, the nor- 
malization of the condensate wave function is defined by 
J dr|$(r)| 2 = N — N, where TV is the total number of 
noncondensate particles. 

As an illustration of this procedure, we consider an 
isotropic harmonic potential V ex t(r) — imw^r 2 , for 
which the order parameter ^(r) is a spherically symmet- 
ric function. In addition, we use the following parameters 



HI): m( s7 Rb) = 1.44 x 1CT 25 kg, v = uj /2tt = 200 
Hz and an s-wave scattering length of a ~ llOao = 
5.82 x 10~ 9 m. Throughout we express lengths and 
energies in terms of the characteristic oscillator length 
d = (h/mtvo) 1 / 2 — 7.62 x 10~ 7 m and the characteris- 
tic trap energy Tiluq — 1.32 x 10~ 31 J, respectively. A 
convenient dimensionless parameter describing the 
effective strength of the interactions is 7 = Na/d, which 
is proportional to the ratio of the average interaction en- 
ergy gn (where n = N/d 3 ) to the characteristic energy 
level spacing htOQ. In the absence of n(r) in Eqs.(lj) and 
(||), the various properties of the condensate scale 
with this single parameter 7. 

For a spherical trap, the eigenfunctions of hg and the 
quasiparticle excitations can be classified according to 
the angular momentum I. The number of excitations 
generated for a given I is determined by the dimension 
of the eigenfunction expansion and is chosen sufficiently 
large to ensure convergence of the noncondensate density 
as determined by the sum over all modes in Eqs.(^) and 
(||). The noncondensate density profile is shown as a 
sum of contributions from each angular momentum I in 
Fig. 1(a) for N = 2000 atoms. The central core of the 
noncondensate is made up from the set of excitations 
with I = 0. 

The fraction of atoms in the noncondensate at T = is 
shown in Fig. 1 (b) as a function of the number of atoms 
N in the droplet. This fraction does not exceed 0.5% of 
the total at T = for the range of iV considered and the 
parameters given above. Having most of the atoms in 
the condensate is in sharp contrast with the situation in 
superfluid He 4 where 90% of the atoms are outside the 
condensate at T = 0. The monotonic increase of N with 
N is due to the dependence of the coupling between the 
u and v functions in the HFB-Popov equations on the ef- 
fective interaction strength 7. The latter increases with 
iV and therefore leads to an increasing number of excita- 
tions out of the condensate. Since iV is a small fraction 
of the total at T = 0, it is to a very good approximation 
a function of the single parameter 7. The value of iV for 
other interaction strengths and harmonic potentials can 
therefore be obtained from the results in Fig. 1(b) by a 
simple scaling with 7. 

With increasing temperature, the noncondensate den- 
sity increases by virtue of the Bose distribution term in 
Eq. (||) . We note that hi (r) depends weakly on tempera- 
ture and is rapidly dominated by 7x2(1") at elevated tem- 
peratures. In Fig. 2(a), we show our self-consistent re- 
sults for h(r) for N = 2000 atoms for a range of tempera- 
tures below the transition temperature of approximately 
100 nK (see below). The insert to Fig. 2 compares n(r) 
and n c (r) at T = 75nK. The trapped atomic gas has 
a two-component structure [p]|l3|] at elevated tempera- 
tures, with a dense core of condensed atoms sitting on 
top of a diffuse gaseous cloud of excited atoms with a 
long tail |16| . Although the noncondensate density in 
Fig. 2(b) looks small in comparison to the condensate, 
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both have an approximately equal number of atoms at 
this temperature because of the r 2 weighting of the inte- 
grated density. We also see that n(r) develops a peak at 
the edge of the condensate with increasing temperature. 
This is the analogue of the sharp peak found in the semi- 
classical Thomas-Fermi approximation [fl3| in which the 
condensate has a sharp cutoff at some radius Rq- 

In Fig. 3(a) we show the total number of atoms N c in 
the condensate as a function of temperature for a 2000 
atom droplet. It can be seen that N c is falling to zero at 
approximately 100 nK, which is close to the BEC tran- 
sition temperature T c predicted by the semiclassical ap- 
proximation |l7],[l8|]. Although there is no sharp transi- 
tion for a finite system, there is still a characteristic tem- 
perature above which N c is small. One way of defining 
this temperature is to fit the temperature variation of N c 
to the functional form N C (T) = A C (0)[1-(T/T C ) Q ], treat- 
ing T c and a as fitting parameters. This functional form 
provides a reasonable parameterization over a wide inter- 
mediate range of temperatures below T c . The values of 
T c extracted in this way are found to scale approximately 
as TV 1 / 3 , which is the iV-dependence found for a trapped, 
noninteracting Bose gas. As can be seen from Fig. 3(b), 
the numerical values of T c found here are slightly smaller 
than the ideal gas values, in a gree ment with recent the- 
oretical [Q and experimental |n| results. The value of 
a extracted from the fits is 2.3, as compared to the ideal 
Bose gas value 18 of 3. 

The size of the noncondensate n(r) is, of course, deter- 
mined by the quasiparticle excitations which deplete the 
condensate. The mode frequencies for the lowest modes 
of angular momentum I = 0, 1 and 2 are shown in Fig. 
4(a) for T = 0. The I = and I = 2 mode frequencies 
split apart from the degenerate harmonic potential eigen- 
value as the effective interaction strength 7 (or N, for 
fixed a/d) is increased. The I = mode is a "breathing- 
type" mode which is influenced by the compressibility 
of the condensate, and increases in frequency as a result 
of the repulsive interactions between the atoms. On the 
other hand, the / = 2 quadrupole mode is a shape reso- 
nance controlled by the surface tension of the drop and 
is seen to decrease in frequency with increasing N. Due 
to the small magnitude of the noncondensate, these fre- 
quencies are very similar to the recent results at T = 
of Singh and Rokhsar Stringari and others 0] in 
which the noncondensate was ignored. 

The lowest lying I = 1 mode is the center-of-mass mode 
of the condensate. It lies very close to, but not exactly 
at, luq. According to t he g eneralized Kohn theorem for 
parabolic confinement P0j| , one would expect a I = 1 
mode at precisely loq corresponding to a rigid oscillation 
of the total droplet density. This property is satisfied 
in the Bogoliubov approximation J5|,^| in which all parti- 
cles are in the condensate. However, in the HFB-Popov 
approximation, the condensate is effectively moving in 
the presence of the 'external' potential V ext (r) + 2gn(r), 
which deviates from the ideal parabolic form. It is only 
because 2gh(v) is a small perturbation in the present case 



that the generalized Kohn theorem is so nearly satisfied. 
An improved approximation in which the dynamics of the 
noncondensate is treated on an equal footing with that 
of the condensate is needed to recover the true center- 
of-mass mode in which both components of the density 
oscillate together. 

We have also examined the quasiparticle excitations as 
a function of temperature. Due to convergence difficul- 
ties at the higher temperatures ]lq ], we show in Fig. 4(b) 
results for a smaller droplet containing N = 200 atoms, 
which has a transition in the middle of the accessible tem- 
perature range. The behaviour for larger drops is qual- 
itatively similar ^3|. For temperatures below the BEC 
transition, the Bose broken symmetry leads quite gener- 
ally to the equivalence of the single-particle and density 
fluctuation spectra However, in the normal phase, 
this equivalence no longer holds and the HFB single- 
particle excitations shown in Fig. 4(b) above T c cannot 
be identified with the density fluctuation spectrum. We 
observe that the splitting between the / = and I = 2 
modes decreases as n c (r) goes to zero through the transi- 
tion temperature, and both modes tend to the harmonic 
excitation frequency of 2wo- The / = 1 dipole excitation 
remains close to u>o for all temperatures. The explanation 
of the high temperature behaviour is straightforward. 
Once n c (r) goes to zero, the single-particle excitations re- 
duce to the normal HF excitations for the self-consistent 
potential V ex t(r) + 2gn(r). Since the total density n(r) 
becomes spread out at temperatures above T c , the HF 
interaction potential provides a fairly constant potential 
energy shift to the harmonic confining potential and the 
lowest excitation energies reduce to those of a noninter- 
acting gas in a harmonic well. 

In conclusion, we have presented a detailed study of the 
behaviour of the condensate and noncondensate density 
profiles as a function of temperature within the Popov 
approximation to the full HFB equations of motion ||. 
Although we have only considered isotropic parabolic 
wells, our method based on Eqs.(|4|)-(||) can also be used 
for anisotropic, as well as anharmonic, traps. Previous 
theoretical studies 13,2^Jl^| have made use of simpli- 
fied semiclassical approximations, whereas our results are 
based on a fully self-consistent numerical treatment of 
the HFB-Popov equations. We have also determined the 
temperature dependence of the excitations through the 
BEC transition and find an interesting evolution from 
collective to single-particle-like behaviour. Our work em- 
phasizes the need for more experimental studies in this 
temperature range. A study of the full HFB equations 
including the off-diagonal correlations neglected in the 
Popov approximation will be reported elsewhere. 

This work was supported by grants from the Natural 
Sciences and Engineering Research Council of Canada. 
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FIG. 2. (a) Noncondensate density for 2000 rubidium 
atoms at various temperatures. (b) The noncondensate 
(solid) and condensate (dashed) densities at T — 75 nK. 



FIG. 3. (a) The number of atoms within the condensate 
for N = 2000, as a function of temperature, (b) The criti- 
cal temperature vs N as calculated self-consistently using the 
HFB-Popov equations (solid line) compared to the ideal Bose 
gas result T c = (haj /k B )(N/((3)) 1/3 (dashed line). 



FIG. 4. (a) Lowest mode frequencies (in units of the trap 
frequency too) for angular momenta I — 0, 1 and 2, as a func- 
tion of the number of atoms in the trap at T = 0. (b) Tem- 
perature variation of the frequencies for N — 200. The lower 
part of the panel shows N c vs T. 



FIG. 1. (a) Noncondensate density at T = for 2000 ru- 
bidium atoms in an isotropic parabolic trap. The innermost 
curve corresponds to the / = excitations, and each successive 
curve includes the contribution of the next angular momen- 
tum. The insert (b) shows the variation of the noncondensate 
fraction (in percent) at T — as a function of N. 
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